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Abstract 

We investigate the approximate dynamics of several dif- 
ferential equations when the solutions are restricted to a 
sparse subset of a given basis. The restriction is enforced 
at every time step by simply applying soft thresholding 
to the coefficients of the basis approximation. By reduc- 
ing or compressing the information needed to represent 
the solution at every step, only the essential dynamics are 
represented. In many cases, there are natural bases de- 
rived from the differential equations which promote spar- 
sity. We find that our method successfully reduces the 
dynamics of convection equations, diffusion equations, 
weak shocks, and vorticity equations with high frequency 
source terms. 



1 Introduction 

In this work, we investigate the approximate dynamics of 
various PDEs whose solutions exhibit behaviors on multi- 
ple spatial scales. These scales may interact with one an- 
other in a non-linear manner as they evolve. Many physi- 
cal equations contain multiscale (as well as multiphysics) 
phenomena, such as the homogenization problems from 
material science and chemistry and multiscale systems in 
biology, computational electrodynamics, fluid dynamics, 
and atmospheric and oceanic sciences. In some cases, the 
physical laws used in the model can range from molecu- 
lar dynamics on the fine scale to classical mechanics on 
the large scale. In other cases, the equations themselves 
contain high wave number oscillations that separate into 
discrete scales, on top of the smooth underlying behavior 
of the system. 

The main source of difficulty in multiscale computation 
is that accurate simulation of the system requires all phe- 



nomena to be fully resolved. The smaller spatial scales 
influence the global solutions, thus they cannot be ignored 
in the numerical computation. In some cases it is possi- 
ble to derive an analytic equation for the effect of small 
scales on the solution (see EQ3]]). In practice, how- 
ever, it may not be possible to derive a simple expres- 
sion that represents the fine scale behavior. Many problem 
dependent methods have been proposed in the literature, 
while a few provide a general methodology for model- 
ing the macroscopic and microscopic processes that yield 
multiscale models. For example, some general meth- 
ods include the heterogenous multiscale method Q, the 
equation- free method ifTUl . multiscale methods for ellip- 
tic problems lTl2lL and the sparse transform method EJ. 
For an overview of general multiscale approaches see []6] . 
A key difference between our method and the methods 
in 1171110111211 is that we are directly resolving all of the sig- 
nificant scales in the solution. By contrast, the methods 
of (71[T0l[T2) directly resolve only the coarse scales of the 
solution, and they separately "reconstruct" the fine scale 
solution (as well as its effect on the coarse scales). 

From the perspective of mathematics, multiscale meth- 
ods began with representation of a function using a global 
basis, such as Taylor series or Fourier series. More so- 
phisticated bases have appeared; for example, any one 
of the many wavelet bases used in imaging and compu- 
tational physics. The key to the basis approximation is 
that each basis element represents behavior on a specific 
scale, therefore the coefficients of the basis provide com- 
plete information about the underlying function. This is 
also the principle behind multiresolution and decomposi- 
tion methods. 

As the methods of multiscale and multiphysics mod- 
eling developed over the past few decades, so did cor- 
responding methods in imaging and information science. 
One of the fundamental ideas in imaging is that of spar- 
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sity. Sparse data representation is used throughout imag- 
ing from compression to reconstruction. Early advances 
in sparse techniques were, e.g., in |[3][5), which presented 
a convex minimization approach to the computationally 
challenging sparse basis pursuit problem. Many mod- 
els have been proposed which use sparsity to produce 
both more efficient numerical methods and better qual- 
ity solutions. Some applications of sparsity to imaging 
include: compressive sensing, reconstruction of images 
from sparse data (EJUS) , and recovery of images using 
sparse regularization ||9][14). The underlying principle of 
sparsity is that images can be approximately represented 
by a small number of terms with respect to some basis. 
Inducing sparsity, creating effective bases, and developing 
efficient computational algorithms have been intensely ac- 
tive fields in information science. 

For imaging and information science, one of the rea- 
sons for the success of sparse methods is their ability 
to resolve drastically different phenomena with a small 
amount of information. This is also a principal goal of 
multiscale modeling. In this work, we transfer sparsity 
methodology, which was developed for information sci- 
ence, to multiscale nonlinear differential equations and 
show that it can be an effective tool for accurately com- 
puting solutions using less information. 

In particular, we propose solving PDEs with the con- 
straint that the approximate solution resides on a sparse 
subspace of a basis. In this way, the complexity of the 
method will depend on the number of basis terms retained 
and (nearly) independent of the grid size. In the following 
sections we will discuss the general problem and the op- 
timization method used to induce sparsity in the solution. 
Also, the general numerical method will be explained as 
well as results of numerical experiments. The method is 
tested on an advection equation with oscillatory velocity, 
a parabolic equation with oscillatory coefficients, a con- 
servation law with oscillatory diffusion, and the vorticity 
equations with high frequency source terms. We conclude 
with a discussion on the proposed work and implications. 

2 Sparsity 

2.1 Problem Statement 

In general, the problem can be stated as follows. Assum- 
ing x G M n and t > 0, let u(x,t) : ft H> K be the 
approximate solution of 

(.» 

subject to the constraint: u(x,t) = ^2jUj(t)(j>j(x), 
where the number of non-zero Cj at a given time step is 



sparse. The operator F(-) can be non-linear, non-local, 
and dependent on the derivative of u. The basis terms <j>j 
are assumed to exist on separate scales, which is true of 
most bases (for example, Legendre polynomials, Fourier, 
Wavelets, etc.). In this way, the basis terms represent dif- 
ferent global behaviors. The method involves two steps: 
evolve the PDE forward in time and project the updated 
solution onto a sparse subset. We first address sparsity 
induction through soft thresholding. 

2.2 Sparsity via Optimization: Soft Thresh- 
olding 

At a given time step, the problem of projecting the up- 
dated solution onto a sparse subset is equivalent to fitting 
a solution u n with corresponding coefficients {u 7 -} at the 
nth time step to a solution u whose corresponding coef- 
ficients {uj} are sparse. This can be written as a con- 
strained least squares fit as follows. 



min \\u - u n \\ 2 L 2 s.t. 



(2) 



& {uj} is sparse 



Expanding with respect to the basis and assuming that the 
basis is orthonormal, this constrained optimization prob- 
lem is related to the following unconstrained problem: 



(L0) 



u = argmin A||w||oH — ||w- 

i7. ^ 



(3) 



where u is the vector of coefficients. The "norm" || • ||o 
is the number of non-zero coefficients in equation [2] This 
makes equation [3] both non-convex and difficult to solve. 
By replacing the L° norm by the L 1 norm, we get the 
following convex relaxation of equation [3j 



(LI) u = argmin A||S||i H — — S n ||^ 2 
it. £ 



(4) 



Note that since u G C, the L 1 norm is H^Hl 1 := J2j 



where \uj\ = y Uj Uj . The solution of equation|4]is given 
by the following equation: 



— max (\u r - 



A,0 



a? 



(5) 



In general, this can be computed for a non-orthonormal 
basis, which is equivalent to a basis pursuit problem with 
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the L 1 norm as a sparse regularizes In that case, the so- 
lution must be found by an iterative method rather than 
the simple shrinkage provided here as an example. The 
resulting minimizer u is a proximal solution which lies on 
a sparse subset of the original coefficient domain [2 ]. This 
can be used to show that the solutions form a contraction 
map in the L 2 norm. Alternatively, we can simply apply 
the soft thresholding on the coefficients directly in order 
to induce sparsity in this way. 

3 Numerical Method 

Assuming u(x, t) is periodic in the domain Q C M n , one 
natural basis is the Fourier basis, whose coefficients are 
the Fourier transform of u{x,t). This is appropriate for 
the examples shown here. For the rest of this work, we 
will use the Fourier basis; however, the overall method- 
ology presented here is independent of the corresponding 
basis. 

Taking the Fourier transform of the PDE from equa- 
tion [T] and discretizing the resulting differential equation 
in time yields a multistep scheme. Since our method does 
not depend on the choice of numerical updating, we can 
assume that the scheme takes the following form: 

v = Q(u n -\...,u n ) (6) 

The updated solution v may be sparse depending on both 
the PDE and the update operator Q, but in general will 
have non-trivial values everywhere depending on the ap- 
proximation and implementation. The auxiliary variable 
v is projected onto a sparse subspace by the shrinkage op- 
erator: 

u n+1 = «S A (v) (7) 
Altogether, the update in the spatial domain is simply: 

U n+1 = J2S X (Q {v n ~\...,v n )) <f>j (8) 

3 

Unlike traditional projections, this is non-linear and adap- 
tive. Rather than sorting the coefficients and retaining a 
fixed number of large amplitude terms or keeping terms 
whose wavenumbers are below some cutoff, the shrink- 
age allows the number and choice of non-zero coefficients 
to evolve over time. Also, this is not the same as hard 
thresholding the solution at every step {i.e. keeping only 
the terms larger than a fixed value) since the coefficients 
that remain have decreased their magnitude by A. Most 
importantly, the projection does not favor any particular 
part of the spectrum; instead the amplitude of the coef- 
ficient determines if it remains. In terms of the Fourier 



basis, the importance is placed on the amplitude not the 
wavenumber. 

For general convergence, as long as A = Cdt p for p 
larger than the accuracy of the scheme used to update the 
variable in time, then the shrinkage operation does not 
change the spatial accuracy of the original method and the 
method will still converge as dt — ^ 0. For example, dis- 
cretize using the forward Euler method, and then expand 
the shrinkage operator to get 

u n+1 = u n + dt F(vF) + O (A) . (9) 

Therefore, to have convergence as dt — » 0, the shrinkage 
parameter must be A = Cdt 1+a . In general, the shrink- 
age operator is non-expansive in each coefficient, hence 
non-expansive in coefficient norms. This may help with 
obtaining a general convergence result. 

4 Numerical Results 

In this section, we discuss the application of the proposed 
sparse method to several equations with different numeri- 
cal schemes. 

4.1 Convection 

The convection equation we consider is the following: 

d t u = a(x)d x u (10) 

where the coefficient a(x) is highly oscillatory. 

Let k be the wavenumber and use spectral Leap Frog as 
the updating for equation [TO1 to obtain 

£n+l = gn-l + 2dt a*(iku U ) (11) 

in which * is the convolution operator over frequency. 

The time step is 0{dx) to preserve the stability condi- 
tion in equation [TO) In Figure [T] the coefficient is chosen 
as follows. 

, N 1 /0.6 + 0.2cos(x)\ 

a(x) = - exp - . — . 

v ; 4 F V 1 + 0.7sin(64x) / 

This choice of a(x) exhibits both fast and slow modes, but 
the particular structure is not directly needed. 

Figure [T] illustrates the performance of the sparse solu- 
tion method on this example by comparison of the sparse 
solution, the true solution produced using a standard fully 
resolved method, and a "low frequency solution" pro- 
duced by solving equation [TT] for wavenumbers k in the 
interval \k\ < K in which K is the number of modes in 
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(a) True (black) and Sparse (blue) Solution in x space (b) True (black) and Sparse (blue 'x') Solution in x space, 

zoomed in 




(c) Sparse (blue) Solution and low frequency (red) Solu- (d) True (black) Solution and Sparse (blue) Solution in u 

tion in x space domain 

Figure 1: Convection with highly oscillatory coefficients. The solution is shown in x space, zoomed in version in x 
space, and in the u space. The solutions are shown on a 512 grid, with dt = 2e — 3, dx = 1.23e — 2, and A = 5e — 03. 



sparse solution. In Figure [T](a-b) the sparse solution pro- 
duced by our method and the true solution at t = 1 are 
plotted in the spatial domain at a given time. In Figure [T] 
(c) the sparse and low frequency solutions are displayed. 
The low frequency solution contains parasitic modes com- 
mon to the Leap frog scheme. In Figure Q](d) the sparse 
and true spectra are plotted. The sparse spectrum captures 
the largest amplitude coefficients throughout the domain. 
In fact, out of the 512 coefficients used in the true solu- 
tion, only 27 are retained in the sparse one (about 5.3%). 

4.2 Parabolic 

The parabolic equation we consider is the following: 

d t u = d x (a(x)d x u) (12) 

where the diffusion coefficient a(x) is highly oscillatory. 
The coefficient is assumed to be bounded, i.e. Am > 



a(x) > A m > 0. This is also related to the elliptic 
case d x (a (x) d x u) = /, since an elliptic equation can 
be solved by taking a parabolic scheme to steady state. 
Alternatively, the corresponding parabolic scheme can be 
iterated forward for a small number of time steps in or- 
der to find a partial solution to the elliptic problem. Then 
by using the partial solution, the locations of the non-zero 
coefficients can be extracted and the elliptic problem can 
be solved by a Galerkin method on these coefficients (4). 

The updated scheme we use for equation [12] is forward 
Euler. 



=u n + dtik a*(iku n ) (13) 

The time step is 0(dx 2 ) to preserve the stability con- 
dition as well as the highly oscillatory nature of the coef- 
ficient a(x) in equation [121 In Figure [2 the coefficient is 
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(a) True (black) and Sparse (blue) Solution in x space 
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(c) True (black) Solution and Sparse (blue) Solution in u 
domain 



(b) True (black) and Sparse (blue 'x') Solution in x space, 
zoomed in 




(d) Diffusion coefficient in x space 



Figure 2: Parabolic diffusion with highly oscillatory coefficients. The solution is shown in x space, zoomed in version 
in x space, and in the u space. The solutions are shown on a 2048 grid, with dt — 1.5e — 8, dx — 3.1e — 3, and 
A = 2.5e-06. 



chosen as follows. 

1 / 0.6 + 0.2cos(» \ 
a(,) = -exp^ + o7sin(256x) j 

In Figure [2] (d) the highly oscillatory diffusion coeffi- 
cient a(x) is plotted in space. In Figure [2] (a-b) the sparse 
solution produced by our method and the true solution at 
t = 1 are plotted in the spatial domain at a given time 
and are nearly indistinguishable. The high frequency in- 
formation is near the scale of the grid size, which can be 
seen in the zoomed in plot. In Figure [2] (c) the true and 
sparse spectra are displayed. The sparse spectrum cap- 
tures the largest coefficients throughout the domain and 
not just the low wavenumbers. In fact, out of the 2048 
coefficients used in the true solution, only 53 are retained 
in the sparse one (about 2.6%). In time, this number of 
non-zero coefficients as well as the identities of the non- 
zero coefficients will change in order to capture various 
behaviors. 



4.3 Viscous Burgers 

To investigate the sparse dynamics of conservative laws 
with diffusion, we use the viscous Burgers type equation. 

d t u + -d x (u 2 ) = d x (a(x)d x u) (14) 

The LHS of equation[l4]is the standard Burgers advection 
term and the RHS is diffusion related to equation [12] The 
equation exhibits three separate phenomenon: 1. Smooth 
large scale behavior from the diffusion, 2. Small scale 
oscillations induced from the high frequencies in the co- 
efficient a(x), and 3. Non-linear interactions between fre- 
quencies from the advection term. The update scheme in 
time is the standard total variational diminishing Runge- 
Kutta 2. 

U! = u n + dtik(a*(iku n )- FfV^) 
v n+1 = i(?T + £i)+ * ifc(a*(ifcSi)-F(^)) 
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where F(u) = \u 1 . As before, we have the stability con- 
dition dt is 0(dx 2 ). 

For Figure [3] the diffusion coefficient is chosen as: 



a(x) = 0.075 exp 



0.65 + 0.2 cos(» 
l + 0.7sin(128x) 



The convolutions in the diffusion and non-linear terms are 
done in the spectral domain, rather than by other meth- 
ods such as the psuedo-spectral method. The various dy- 
namics can be seen in the spatial and in the spectral plots 
(see Figure 0. The true, sparse, and low frequency solu- 
tions are plotted in the space in Figure Oa-b). The low 
frequency projection is done by thresholding any coeffi- 
cients outside of a particular range. Specifically, the num- 
ber of low wavenumbers retained is the same as the sparse 
solution, although their identities are dramatically differ- 
ent. The sparse solution captures the local and global be- 
haviors of the solution more accurately than the low fre- 
quency projected solution. In Figure [3] (c-d), the spectrum 
of the true solution is compared to the sparse and low fre- 
quency spectra, respectively. The local peaks in the spec- 
tra are related to the wavenumbers in the diffusion coef- 
ficient a(x) and the harmonics induced by the non-linear 
advection term. Notice that in this case, each of the dis- 
tributions in the spectral domain do not decay as rapidly 
as in the parabolic case. The sparse solution contains 130 
coefficients out of a total possible 1024, about 12.7%. 

4.4 Vorticity Equations 

The vorticity equation we consider is derived from the (2- 
d) incompressible Navier-Stokes equation [11]: 



d t u + (\/ ± A~ 1 u) • \7u = <yAu + / 



(15) 



where u is the vorticity (not the velocity). Similarly 
to equation [H equation [15] exhibits three separate phe- 
nomenon: 1. Smoothness from the diffusion term on the 
RHS, 2. Small scale oscillations induced from the high 
frequencies in the source term /, and 3. Non-linear inter- 
actions between frequencies from the advection term on 
the LHS. However, since the operator V x A -1 is smooth- 
ing (in some sense), the advection term can be viewed as 
less non-linear than the one found in equation[l4] In terms 
of the numerical method, the operator V^A -1 dampens 
the coefficients by a factor which acts as |fc| _1 . 

For the numerical implementation, the diffusion term 
is discretized using Crank-Nicolson while the advection 
term is lagged. Since the operators are diagonalized in the 
coefficient basis, the steps can be invertible and lead to a 
simple updating scheme. 



v n+1 



2dt 



2 + 7<ft|fc| 2 
2-jdt\k\ 2 



(ik 1 - (\k\~ 



* iku 11 



■n 



2 + 1 dt\k\ 2 
For Figure HI the forcing term is chosen to be: 



f(x,y) = 0.025 



sin(32x) + sin(32y) 
1 + 0.25 (cos(64r) + cos(64y)) 



The standard stability condition is used for choosing the 
time steps in order to insure capture of all small scale be- 
haviors. In Figure [4] (a-b), the true and sparse solutions 
are plotted in the spatial domain. Notice the oscillations 
introduced by the source term interact with the two vor- 
tex patches and thus contribute to the global behavior of 
the solution. The spectra of the true and sparse solution 
is plotted in Figure [4] (c-d), where the low wavenumbers 
are located in the middle of the domain. The sparse solu- 
tion retains about 3.95% of the coefficients. In the sparse 
spectrum, the coefficients are located throughout the do- 
main, including the highest frequency itself (seen on the 
boundary of the spectral domain). In Figure [4] (e), the L 2 
and L°° error are shown to decrease as the resolution in- 
creases. This sparse solutions, as well as the other exam- 
ples presented here, converge as the spatial discretization 
goes to zero. 

It was observed that as the dimension increases, the 
sparsity of the solution also increases (since it is propor- 
tional to the product of the sparsities in each dimension). 
Thus the method scales well with dimension. 

It is worth noting that in a related work, wavelet hard 
thresholding was used to separate coherent and incoherent 
structures in turbulent flows JSJ. 



5 Discussion 

In the examples from the previous section, the PDEs con- 
tained a mixture of multiscale properties with a diffusion 
term. The combination of non-linear and oscillatory terms 
created a large range of wavenumbers in the solution, 
while the dynamics produced a range of amplitudes. This 
gives the necessary structure for sparsity with respect to 
the Fourier basis. In general, the highest order derivative 
will determine the appropriate basis in which the solutions 
could be sparse. 

If the spectrum is more localized, i.e. non-zero regions 
in the low frequency regime, then the proposed model can 
better condition the numerical method. Empirically, the 
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(a) True (black), Sparse (blue), and Low Frequency (b) True (black), Sparse (blue), and Low Frequency 

(red) Solution in x space (red) Solution in x space, zoomed in 




(c) True (black) and Sparse (blue) solutions in u do- (d) True (black) and Low Frequency (red) solutions in 

main u domain 

Figure 3: Viscous Burgers equation with highly oscillatory diffusion. The solution is shown in x space, zoomed in 
version in x space, and in the u space. The solutions are shown on a 1024 grid, with dt = 7.6e — 6, dx = 6.1e — 3, 
and A = 6.3e — 5. 



shrinkage operator acts as a non-linear filter on the coef- 
ficients. It was observed that for a fixed C and p, where 
A = Cdt p , the numerical updates presented here with dt 
larger than theoretically and numerically possible in the 
original scheme will converge. In the case of the vorticity 
example, dt can be taken much larger when soft thresh- 
olding is applied than in the standard scheme. Also, the 
non-linear filter seems to reduce parasitic modes and spu- 
rious oscillations found in spectral approximations for lin- 
ear and for non-linear slightly viscous hyperbolic equa- 
tions (see Figures Q](b) and (c)) and 0(g). 

One key point is that our method works by fully re- 
solving the solution. Its efficiency is gained by omitting 
modes that are insignificant. This requires that A is small 
enough that the filter does not annihilate essential fea- 
tures. For example, if the initial data is smaller than A 
for a particular unstable mode, then our approximate so- 
lution will not match the true dynamics. As the grid is 
refined, the mode will be captured (since A decreases as 
Ax decreases). 

In terms of complexity, each iteration is dominated 



by the convolution step. The convolution in the coef- 
ficient domain (spectral domain) can be done explicitly 
over the n s (t) -sparse vectors rather than transforming 
back onto the spatial grid, which is O (n s (£) 2 ) at each 
step. When n s (t) 2 << NlogN, convolving in the spec- 
tral domain rather than transforming back and forth be- 
tween domains decreases the computational cost of the 
algorithm. Knowing a priori the maximum sparsity, i.e. 
n s ,max = max t n s (t), faster routines and transforms 
could be optimized for specific problems and applications. 
For example, one can optimize the transform between the 
spatial and coefficient domains knowing the given sparsity 
at the current step and the non-trivial coefficients' identi- 
ties. In the linear cases, as in equation Q31 the operation 
a* can be stored as a large but sparse matrix, reducing 
the updates to a sparse matrix - sparse vector operation 
at every iteration. Our goal in this work is to formulate 
a PDE solver that promotes sparsity. In future work we 
will present a study of the computational complexity and 
speed. 

When the dynamics are dominated by a linear term, for 
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(a) True Solution in x space (b) Sparse Solution in x space (c) Low Frequency Solution in x space 




step size 

Figure 4: Vorticity equations with high frequency source term. The solution is shown in x space, zoomed in version in 
x space, and in the u space. The solutions are shown on a 256 by 256 grid, with dt = 0.025, dx = 0.0245, 7 = 0.001, 
and A = 0.0497. 



example high viscosity, the identities of the non-trivial 
coefficients settle over time. This was also observed in 
the non-linear cases, but over a longer time period. This 
would enable creation of a sparse basis for elliptic equa- 
tions (e.g., with oscillatory coefficients). 

In many of the cases here, it is possible to get hyper- 
sparse solutions (those with 1 % or fewer coefficients) at 
the cost of accuracy. 

6 Conclusion 

In this work, we have proposed a method to fully resolve 
the solutions of multiscale PDEs while only retaining im- 
portant modes. The reduced dynamics created by the 
sparse projection properly captures the true phenomena 
exhibited by the solution. The sparse projection amounts 
to a shrinkage of the coefficients of the updated solution at 
every time step. There exist many possibilities for using 
the sparsity to optimize individual algorithms and create 



faster, more efficient computational routines. 
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